%% First load data arrays for analysis
fnames = dir("monod*.mat");
for k=1:length(fnames)
    load(fnames(k).name);
end

%% Define molar mass
clc;
close all;
%% Do 25C
tempmat = monod_mat_proline_25C;
smat = size(tempmat);
ind = 1:3;

conc1 = tempmat(ind,1); % Convert to mM
grmax1 = nanmean(tempmat(ind,2:3),2);
grmax_err1 = nanstd(tempmat(ind,2:3),0,2);

% Do fit on corrected concentration (in mM)
figure;
errorbar(conc1,grmax1,grmax_err1, 'bo', 'linewidth', 1);
    
hold on;

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')
box off;

% Monod fit
[f1,f2] = monod_fit_weighted(conc1,grmax1,grmax_err1);
cplot = linspace(0,max(conc1),100);
plot(cplot,f1(cplot),'b');
ci = confint(f1);

KM = f1.b;
KMerr = abs(ci(2,2) - KM);

disp('Fit and error of KM (mM)');
disp(KM);
disp(KMerr);


%% Do 30C
tempmat = monod_mat_proline_30C;
smat = size(tempmat);

conc1 = tempmat(:,1); % Convert to mM
grmax1 = nanmean(tempmat(:,2:3),2);
grmax_err1 = nanstd(tempmat(:,2:3),0,2);

% Do fit on corrected concentration (in mM)
hold on;
errorbar(conc1,grmax1,grmax_err1, 'ko', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')
box off;

% Monod fit
% Monod fit
[f1,f2] = monod_fit_weighted(conc1,grmax1,grmax_err1);
cplot = linspace(0,max(conc1),100);
plot(cplot,f1(cplot),'k');
ci = confint(f1);

KM = f1.b;
KMerr = abs(ci(2,2) - KM);

disp('Fit and error of KM (mM)');
disp(KM);
disp(KMerr);

%% Do 37C
tempmat = monod_mat_proline_37C;
smat = size(tempmat);

conc1 = tempmat(2:end,1); % Convert to mM
grmax1 = nanmean(tempmat(2:end,2:3),2);
grmax_err1 = nanstd(tempmat(2:end,2:3),0,2);

% Do fit on corrected concentration (in mM)
hold on;
errorbar(conc1,grmax1,grmax_err1, 'ro', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
[f1,f2] = monod_fit_weighted(conc1,grmax1,grmax_err1);
cplot = linspace(0,max(conc1),100);
plot(cplot,f1(cplot),'r');
ci = confint(f1);

KM = f1.b;
KMerr = abs(ci(2,2) - KM);

disp('Fit and error of KM (mM)');
disp(KM);
disp(KMerr);